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ABSTRACT 

We apply our recently suggested statistical approach to thermal evolution of isolated neu¬ 
tron stars and accreting quasi stationary neutron stars in X-ray transients for constraining the 
position and relative broadening a of the direct Urea threshold of powerful neutrino emis¬ 
sion in neutron star cores. We show that most likely explanation of observations corresponds 
io a ^ 0.08 - 0.10 and to the neutron star mass, at which the direct Urea process is open, 
Md ~ 1.6- 1.8 Mo. 
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1 INTRODUCTION 

It is well known that theories of thermal evolution of neutron 
stars have potential to explore properties of superdense matter in 
neutron star cores ( [Page et ^|2004| [2009 1 [Yakovlev & Pethic^ 
2004 jl. Recently w e have suggested a statistical approach ( jBezno- 
gov & Yakovlev||2015[ hereafter Paper I) to study thermal evo¬ 
lution of middle aged {t ~ 10^-10^ yr) isolated neutron stars 
(INSs) and old {t > 10^ yr) transiently accreting quasistation¬ 
ary neutron stars in low-mass X-ray binaries (in X-ray transients; 
XRTs). INSs cool down loosing their thermal energy; neutron stars 
in XRTs are warmed up by deep crustal he ating (jHaensel & Zdunikj 
|1990| [2008} [Brown, Bildsten & Rutledge||1998| ) during accretion 

episodes. Thermal evolution of INSs and neutron stars in XRTs 
allows one to test the same physics of superdense matter in neu- 
tron star cores (e.g., [Yakovlev & Haensel|2003] [Yakovlev, Leven-j 

[fish & Haensel|2003| ). Traditional theories of thermal evolution of 

these objects are based on calculations of cooling curves for INSs 
(e.g., T^{t)), and heating curves for accreting neutron stars (e.g., 
L“((M))). Here, is the redshifted effective surface temperature 
of the star, t is the stellar age, is the redshifted photon surface 
luminosity of a neutron star in quiescent states of XRT, and {M) is 
the mass accretion rate in XRT averaged over neutron star cooling 
time-scales (about 1 kyr or more). 

The statistical approach of Paper I replaces theoretical cool¬ 
ing/heating curves by probabilities to find neutron stars in certain 
regions of the - t or - (M) planes. These probabilities are 
obtained by averaging individual cooling/heating curves over mass 
distributions of isolated or accreting neutron stars and over distribu¬ 
tions of other neutron star parameters, particularly, of mass of light 
elements in the heat blanketing envelopes. Light elements increase 
heat transparency of the heat blanketing envelopes in comparison 
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with standa rd iron envelopes (e.g.,[Gu dmun dsson, Pethick & Ep-| 
[stein[[1983[|Potekhin, Chabrier & Yakovlevnl997[ l and affect thus 
the thermal evolution of neutron stars. Taken an increasing statis¬ 
tics of the sources, calculated probabilities simplify comparison of 
the theory with observations. Moreover, these probabilities depend 
not only on the properties of superdense matter in neutron stars but 
also on distributions of neutron star parameters (first of all, on the 
mass distributions) which enables one to study these distributions 
together with the properties of superdense matter (Paper I). 

Full implementation of the statistical approach is a compli¬ 
cated task. In Paper I the idea was illustrated using neutron star 
models with nucleon cores having one equation of state (EOS) of 
superdense matter. This EOS was constructed by [Kaminker et al.[ 
( [2014| ) and abbreviated as HHJ because it belongs to the family of 
EOSs suggested by |Heiselberg & Hjorth-Jens^ ( [1999[ ). A prelimi¬ 
nary analysis of Paper I indicated that observational data on isolated 
and accreting neutron stars could be reconciled if powerful direct 
Urea process of neutrino emission ( [Lattimer et al. 11991] ) operated 
in the cores of massive stars and the threshold for the onset of the 
direct Urea process with growing density p were broadened (e.g. 
[Yakovlev et al.|200T) . Otherwise the theory cannot explain a num¬ 
ber of observed sources. Typical masses of INSs need to be lower 
than the mass Mq at which the direct Urea process is on, but accret¬ 
ing neutron stars should be overall more massive to explain at least 
one source, SAX Jl808.4-3658 ([Campana et al.[2002|[Galloway &] 
[Cumming[2006[[Heinke et al.[2009[ >, where the direct Urea process 

seems to operate. 

This paper extends the analysis of Paper I in two ways. Eirst, 
we consider two EOSs of superdense matter instead of one which 
illustrates the effects of EOSs. Secondly, we take into account that 
one should distinguish a formal threshold of the direct Urea pro¬ 
cess (at a density p - poo at which the process becomes allowed 
by momentum conservation of reacting particles for a particular 
EOS; [Lattimer et al.|1991[ l and the actual threshold Pd which can 
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Table 1. Middle-aged cooling INSs (left) and accreting neutron stars in 
XRTs (right) whose thermal surface emission has been detected or con¬ 
strained; see Paper I for details. 


Table 2. Gravitational masses M, central densities pcu (in units of 10^^ 
g cm“^) and circumferential radii R of neutron star models with the HHJ 
andBSk21 EOSs. 


Num. 

INS 

Num. 

XRT 

1 

PSRJl 119-6127 

1 

Aql X-1 

2 

RX J0822-4300 

2 

4U 1608-522 


(in Pup A) 

3 

MXB 1659-29 

3 

PSRJ1357-6429 

4 

NGC 6440 X-1 

4 

PSR B0833-45 

5 

RX J1709-2639 


(Vela) 

6 

IGR 00291-^5934 

5 

PSR B1706-44 

7 

Cen X-4 

6 

PSR J0538-F2817 

8 

KS 1731-260 

7 

PSRB2334-F61 

9 

IM 1716-315 

8 

PSRB0656-F14 

10 

4U 1730-22 

9 

PSRB0633-F1748 

11 

4U 2129-^47 


(Geminga) 

12 

Terzan 5 

10 

PSR B1055-52 

13 

SAX J1808.4-3658 

11 

RXJl 856.4-3754 

14 

XTEJ1751-305 

12 

PSR J2043-F2740 

15 

XTE J1814-338 

13 

RX J0720.4-3125 

16 

EXO 1747-214 

14 

PSRJ1741-2054 

17 

Terzan 1 

15 

XMMU J1732-3445 

18 

XTE 2123-058 

16 

Cas A neutron star 

19 

SAXJ1810.8-2609 

17 

PSR J0357-F3205 

20 

IH 1905-^000 


(Morla) 

21 

2S 1803-45 

18 

PSRB0531-F21 (Crab) 

22 

XTE J0929-314 

19 

PSR J0205-F6449 

23 

XTE J1807-294 


(in 3C 58) 

24 

NGC 6440 X-2 


be shifted with respect to poo- The shift can be produced, for in¬ 
stance, by superfluidity of neutron star matter; superfluidity can 
suppress the direct Urea process and increase po (e.g., [Yakovlev] 
jet al. 12001] [Yakovlev & Pethick|2004| l. In addition, the shift can be 
produced by strong magnetic Adds in neutron stars cores ( [Baiko &[ 
[Yakovlev|1999[[^kovlev et al.|2001| ); the flelds modify momen¬ 
tum conservation rules and can decrease po with respect to poo- It 
is also possible that a threshold for enhanced neutrino emission can 
be broadened and shifted by nuclear physics effects (as discussed, 
for instance, by [Schaab et aL|19^[Blaschke, Grigorian & Voskre-j 
[sensky|2004[ and in references therein). Let us stress that it is the 
actual threshold po which regulates thermal evolution of neutron 
stars; this threshold is basically unknown even if we knew the EOS. 
Therefore, we will treat Pd as a free parameter and try to constrain 
its values. Note that the same effects which shift po can broaden the 
direct Urea threshold. Thus, although we consider threshold broad¬ 
ening and shifting as independent phenomena, it may be actually 
not so. In addition, it is important to stress that large shifts and 
strong broadening seem unlikely from theoretical points of view. 

As in Paper I, we will treat the broadening of the direct 
Urea threshold on phenomenological level by introducing the same 
model broadening function characterized by the relative broaden¬ 
ing width a. We will not use any specific physical broadening 
model (superfluidity, magnetic flelds, nuclear physics effects) but 
postpone such studies till future, more detailed consideration. Our 
present analysis will be qualitative but hopefully it will give gen¬ 
eral understanding of the problem. The observational basis will be 
the same as in Paper I (where detailed list of references is given). 
For convenience of the reader Table [T] lists observational sources 
plotted in the figures below. 


EOS / model 


HHJ 



BSk21 


M/MO 

Pcl4 

R(km) 

M/MO 

Pc 14 

R(km) 

Maximum mass 

2.16 

24.5 

10.84 

221 

22.9 

11.04 

Direct Urea onset 

1.72 

10.0 

12.49 

1.59 

8.21 

12.59 


2 QUALITATIVE ANALYSIS 


We adopt two EOSs, (i) the HHJ EOS used in Paper I and (ii) 
the BSk21 EOS ( Go r iely, Charnel & Pearson|2010| [Pearson et al.] 
[2012[ [Potekhin et al.|2013[ l. The maximum-mass and direct-Urca- 
onset neutron star models for these EOSs are presented in Table 
1^ The important parameter here is the formal (determined by a 
given EOS) mass density threshold poo for the onset of the electron 
direct Urea process and the associated minimum mass of the star 
Mdo in which the direct Urea occurs. As seen from Tablep™ = 
1.00 X lO^^ g cm-^ = 1.72 MG and pgg^^i ^ 3 ^ >< iqu 

g cm"^, = 1.59 Mo. Note that similar table 3 of Paper I 

contains two minor typos which do not affect the results of Paper I 
because calculations were made with correct parameters. 

As in Paper I, the neutrino emissivity Q of the electron direct 
Urea process is presented as 2 = Qob, where Qq is the basic direct 
Urea emissivity disregarding sharp (step-like) threshold at p = poo 
( [Lattimer et al.|199r| ; bis a. phenomenological function to shift and 
broaden the threshold. It is chosen as 


b{x,^) = 


0.5 [1 -h erf (v - ^cxj)] 
0 
1 


at - du < V - < du, 

at X - ^(jj < -du, 

atx-^(jj>(jj. 


( 1 ) 


Here, erf(v) is the error function; x - {p - Pdo)/(q^Pdo). ~ 
^PY)lpm being the threshold broadening factor and Apo a char¬ 
acteristic broadening density interval; ^ (with -1 1) specifies 

the position of the actual electron direct Urea threshold, p = po 
(Section [^. In fact, b is determined by two physical quantities, po 
and dr. At ^ = 0 the threshold is not shifted from the ‘basic’ value 
Pd = Pdo prescribed by EOS. If dr ^ 0, the function b becomes 
step-like with the jump at p = po- If ^ ^ 0, the threshold po is 
shifted. 


Pd =Pdo(1(2) 

where co ^ 5.40988 is an auxiliary quantity which defines the cut¬ 
off level 10 "^"^ of the broadening function, 0.5 [1 - erf(dD)] = 10 “^"^. 
At p = Pd the broadening function b - 0 . 5 ; the function varies 
from 0.1 to 0.9 in the density interval - 0.91 a ^ (p - Pd)/Pd ^ 
0.91 dr. The shape of the broadening function is actually the same 
as in Paper I but the threshold position is allowed to be shifted from 
Pdo- The muon direct Urea threshold is shifted and broadened in 
the same way as the electron one which seems rather unimportant. 
Note that the shape of our broadening factor b as a function of p 
is symmetric with respect to p = pd (fig. 8 of Paper I). In reality, 
it can be asymmetric and the asymmetry may affect evolution of 
neutron stars which remains to be studied in the future. 

Let us perform a qualitative analysis of the effects of the direct 
Urea process in neutron stars. These effects are characterized by 
two parameters, Pd and a. In addition, we introduce minimum and 
maximum masses of neutron stars involved in our analysis, 
and They may be not the absolute minimum and maximum 
masses of neutron stars (e.g., [Haensel, Potekhin & Yakovlev|2007[ l 
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Figure 1. Ranges of possible values of Md versus a (shaded area) for the 
HHJ EOS at = 1.1 Mo and =2.1 Mo. The dashed line shows 
^Dmax while the solid line is MDmin- The double dot-dashed line is 
(Table[g. See text for details. 


but rather minimum and maximum masses in the ensemble of INSs 
and accreting neutron stars in XRTs we would like to study. The 
masses of neutron stars in XRTs should be overall higher because 
of accretion. Therefore, refers to XRTs. 

By varying pd and a, we can estimate upper and lower values 
of Pd (pomax and pomin) and associated upper and lower masses 
(^Dmax and Momin) at which the direct Urea process appears. If 
Pd is too low, the direct Urea is allowed in low-mass stars making 
them cooler. If, on the other hand, Pd is too high, the direct Urea is 
suppressed in high-mass stars making them warmer. 

Both values, MDmax and MDmin, can be estimated by com¬ 
paring theoretical cooling/heating curves with observational data. 
They depend on the broadening factor a and on assumed val¬ 
ues of M^.j^ and M^^^. Our formal procedure to estimate possible 
ranges of a and Pd (or Md) is like this. First, too small broaden¬ 
ing {a < 0.05) does not allow us to explain observations of many 
isolated and accreting neutron stars (fig 9-12 of Paper I). With this 
in mind we consider the case of cr > 0.05. For each a we calculate 
warmest cooling and heating curves (taking M = M^.^ ) at differ¬ 
ent values of pd- The minimum value pDmin is determined as the 
one which shifts the warmest cooling or heating curve (the values 
of or L“) down by 5% with respect to the curve calculated for 
non-broadened direct Urea threshold. Of course, the assumed 5% 
level is conditional and can be changed but with qualitatively the 
same results. It is important that this procedure gives nearly the 
same pDmin (or MDmin) for INSs and XRTs. 

As mentioned above, the coldest neutron stars (M = M^^^) are 
thought to be in XRTs. Their heating curves can also be calculated 
for different values of pd at each a ^ 0.05. We estimate pDmax (or 
MDmax) as such which shifts the coldest heating curve up by 5%. 
This estimate is almost insensitive to a. 

The results of this analysis are presented in Figs. [T]-0Fig.[T] 
shows allowable values of Md and a (shaded domain) for the HHJ 
EOS assuming M^^^ = 1.1 MO and = 2.1 MO at cr > 0.05. 


Figure 2. Allowable minimum (short-dashed and solid lines) and maxi¬ 
mum (long-dashed and dash-dotted lines) masses Md of neutron stars ver¬ 
sus a for the BSk21 EOS at M*. =1.0 and 1.1 Mo and M* =2.1 and 

min iiiax 

2.2 Mo. The entire shaded area shows the range of possible values Md for 
M^in = 1.0 Mo and M^^^ - 2.2 Mo, while the denser shading is the same 
for M^.^ = 1.1 Mo and MJ^^^ = 2.1 Mo. See text for details. 



a 

Figure 3. Comparison of minimum and maximum allowable masses Md 
versus a for the HHJ and BSk21 EOSs at M^.^ = 1.1 Mo and M^^^ = 
2.1 Mo. See text for details. 


The solid line is MDmin and the dashed line is MDmax- The double 
dot-dashed line is Mdo (Table [T}. 

Fig. presents MDmin and MDmax as a function of a for 
the BSk21 EOS. The solid and short-dashed lines are MDmin for 
M^n = 1.0 and 1.1 Mo, respectively. The dot-dashed and long- 
dashed lines show MDmax for M^^^ - and 2.2MO, respec¬ 
tively. Light-greyed and dark-greyed areas (combined together) 
show the wide range of allowed values of Md corresponding to 
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the mass range of 1.0-2.2 MG, while dark-greyed area (alone) 

shows narrower range of allowed Mq values corresponding to the 
1.1-2.1 Mo mass range. 

Fig. [^compares minimum (solid and long-dashed curves) and 
maximum (dotted and dot-dashed curves) allowable masses 
versus a for the HHJ and BSk21 EOSs at = 1.1 Mo and 

^;ax = 2.1MO. 

If Momin formally exceeds Momax. the broadening of the direct 
Urea threshold becomes too strong and theoretical cooling/heating 
curves cannot explain all observational data. The presented analysis 
gives a tighter upper constraint on a than in Paper I. As seen from 
Figs. □-HI a should not exceed 0.12 —0.15. 

Also notice that the derived constraints on Md depend on as¬ 
sumed values of and (which restrict masses of neu¬ 
tron stars of our interest). We have taken = 1.0,1.1 Mo and 
M^ax = 2.1,2.2 Mo which do not contradict observational data. In 
the next section we use = 1.1 Mo and = 2.1 Mo to be 
consistent with Paper I. 


3 POSITION AND BROADENING OF DIRECT URCA 
THRESHOLD 

3.1 Illustrative mass functions 


Let us supplement the qualitative analysis of the previous section by 
some numerical results. We will mainly present probabilities to find 
isolated and accreting neutron stars in certain regions of 7“ -1 (for 
INSs) or - (M) (for XRTs) planes compared with observations 
(Section Table see also Paper I, tables 1 and 2). Specifically, 
they are the differential probabilities 


dPi = p, (T-, t) drr, dPa = p, (l;, (M» di;, (3) 

to find an isolated (i) or accreting (a) neutron star of given age t or 
mass accretion rate (M) in small surface temperature or luminos¬ 
ity intervals, dP“ or dL“, respectively. By construction, the total 
probabilities f pi dP“ and f p^ dL“ are conserved (independent of 
t and (M)). 

The differential probability densities, piiT^j) and 
p^{L^,{M)), will be plotted by greyscaling (in relative units). 
The denser the scaling, the larger the probability. White regions 
refer to zero or very low probability. A graphical example of the 
probability distribution is discussed in the end of Section [33] 

The probabilities are calculated (Paper I) by averaging the 
families of cooling or heating curves over mass distributions of 
isolated or accreting neutron stars and over masses of light ele¬ 
ments AMie in the heat blanketing envelopes of the stars. For clarity, 
by dashed lines we will present also the initial cooling or heating 
‘reference’ curves for neutron stars with masses M - 1.1, 1.2,..., 
2.1 Mo with iron heat blankets. Higher curves correspond to lower 
M. Curves for low-mass stars (M < Md) often merge in nearly the 
same curve. 

As in Paper I, the distribution function over AMie will be taken 
uniform, and AMie/M is allowed to vary from 0 to its maximum 
value 10“^. Following Paper I we use normal mass distributions 
/(M) for INSs and lognormal for neutron stars in XRTs, 


/i(M) = 


1 1 
M ^/2nc^i 


exp 


2o-f /’ 


/a(M) = 


1 1 

Tr~F^ - 

dYrtMa-^ 


(ln[M/M0]-;Uaf 


(4) 


Here ‘i’ stands for INSs and ‘a’ is for XRTs; cti, pi, and cr^, are 


Table 3. Parameters of model mass distribution functions of INSs (i) 
and accreting neutron stars (a) in XRTs and corresponding values of Md 
used in illustrative examples. See also Fig.[^ 


Model 



P^ 


Md 

1 

1.40 

0.15 

0.47 

0.17 

1.72 

2 

1.35 

0.15 

0.42 

0.18 

1.66 

3 

1.47 

0.10 

0.43 

0.09 

1.65 

4 

1.25 

0.10 

0.30 

0.12 

1.45 


the parameters of these distributions; Ni and Aa are normalization 
factors. We will take four pairs of mass functions 1-4 listed in Table 
Note that cti and pi are expressed in units of Mo, while cTa and p^ 
are dimensionless. We consider these distributions within the mass 
range from M^.j^=l.l MO to MO, with /(M) = 0 outside 

this range. These distributions are plotted in Fig. [^ Their parame¬ 
ters have been chosen ‘by eye’ to get better agreement between the 
theory and observations in each specific case described below. The 
distributions li and la were used in Paper I; other distributions are 
new. Evidently, these distributions are not unique; a detailed study 
of most suitable distributions will be done in future publications. In 
all four cases (1-4) the peaks of /a(M) are shifted to higher M (by 
~ 0.15 Mo) with respect to the peaks of as a natural result 

of accretion. By way of illustration, we will also take the uniform 
mass distribution within the same mass interval from 1.1 to 2.1 MO 
for INSs and XRTs (plotted by the dot-dashed curve ‘0’ in Fig.j^. 


3.2 Very small broadening, cr < 0.05 

In Paper I, we demonstrated the effect of too small broadening, 
with a < 0.05. It splits the populations of INSs and XRTs into 
families of rather warm and cold sources separated by ‘gaps,’ in 
disagreement with observations (figs 9-12 of Paper I). It was shown 
for the HHJ EOS. Now we obtain similar results for the BSk21 
EOS. Note that at cr ^ 0.05 the ‘gap’ is narrow and, in principle, 
could be closed by tuning mass functions. However, to this aim one 
needs sufficiently sharp and narrow mass distributions which seem 
unlikely. Therefore, the conclusion that the broadening of the direct 
Urea threshold is not too small looks solid. At cr > 0.05 the ‘gaps’ 
disappear and the theory can be consistent with the data. 


3.3 Most ‘successful’ broadening, a ^ 0.08 —0.10 


As shown in Paper I (for the HHJ EOS), the case of cr ~ 0.1 is most 
suitable to explain the observations. Let us focus on this important 
case. For a fixed cr ~ 0.1 one should find most appropriate value 
of Md and most appropriate mass functions /i(M) and /a(M) to 
explain the data. 

Fig. [^ shows the probability distributions for INSs (in the 
- t plane, left-hand panel) and XRTs (in the - (M) plane, 
right-hand panel) assuming the HHJ EOS, cr = 0.1, Md = = 

1.72 Mo, and the mass functions li and la (after figs 13 and 14 of 
Paper I). This model explains the observations reasonably well. Ac¬ 
cording to Fig.[^ we can slightly vary Md around Mdo and slightly 
adjust the mass functions. These new explanations will also be sat¬ 
isfactory although the explanation of Paper I is one of the best. 
The source 15 on the left-hand panel (XMMU J1732-3445, Table 
[Klochkov et al.|20T3]l is too hot to be explained by the present 
model. As shown in Klochkov et al.[p015| ) and Paper I, it can be 
interpreted as a neutron star with strong proton superfluidity in the 
core and heat blanketing envelope fully made of light elements. The 
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Figure 4. Illustrative models 1-4 (Tablej^ of mass distribution functions for INSs (i, left-hand panel) and accreting neutron stars in XRTs (a, right-hand 
panel). The dot-dashed line ‘0’ refers to the uniform mass distribution. Vertical dashed lines show values of Md for corresponding models. 




Figure 5. Probability distributions for INSs in the -1 plane (left) and for accreting neutron stars in the - (M) plane (right) compared with observations. 
The HHJ EOS is used; the distributions of neutron star masses are li and la, and the distribution over the mass of light elements in surface layers is uniform. 
Dashed lines show 11 ‘reference’ cooling curves for stars with iron envelopes and masses (from top to bottom) M = 1.1,1.2,..., 2.1 Mo. The direct Urea 
threshold is broadened but not shifted, a = O.l, Md = details. 


neutron star in Pup A (source 2) also seems to require proton super¬ 
fluidity and the heat blanket of light elements but the requirement is 
much less stringent. We do not include the effects of strong proton 
superfluidity here to simplify our analysis. 

For illustration, Fig. [^presents also the probability distribu¬ 
tions for isolated and accreting neutron stars with the HHJ EOS 
but using the uniform mass distribution ‘0’ instead of li and la. 
Such a model was not studied in Paper L In contrast to the previous 
one (Fig.|^, it predicts the existence of cold INSs which have not 
been observed. In addition, it is less successful in explaining the 
existence of many neutron stars in XRTs which are intermediate 


between coldest and hottest ones. The presented example demon¬ 
strates the importance of mass distribution functions f(M) in our 
analysis. 

Now let us consider the BSk21 EOS, whose formal direct Urea 
threshold corresponds to Mdo = 1.59 Mg. As seen from Fig. 
the values a = O.l and Md = ^d? most suitable for ex¬ 

plaining the observations. After several trial calculations we have 
found that the case of cr = 0.09 and Md = 1.66 Mo, where the 
actual direct Urea threshold is shifted to higher densities, is more 
successful. The distributions of isolated and accreting neutron stars 
over respective diagrams for this case are plotted in Fig. The 
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log (m'^, Mo/yr 


Figure 6. Same as in Fig.[^but for the uniform mass distributions. See the text for details. 




log Mo/yr 


Figure 7. Same as in Fig.|^but for the BSk21 EOS, a = 0.09, Md = 1.66 Mo, and neutron star mass distributions 2i and 2a. 


mass distributions 2i and 2a (Fig. for this case are also modi¬ 
fied with respect to li and la to achieve better agreement with the 
data. The modified mass distributions contain slightly less massive 
isolated and accreting neutron stars. Otherwise, theoretical models 
presented in Figs.[^and|^look rather similar. Note that sufficiently 
large masses Mq > 1.5 M© agree with some theoretical expecta¬ 
tions (e.g., IKlahn et al.|2006l l. 

To visualize the theoretical probabilities ([^, in Fig.j^we show 
the calculated differential probability distributions dP/dT^ of INSs 
as functions of for our ‘most successful’ models. Here, is 
the effective surface temperature expressed in MK, so that dP/dT^^ 
is the probability to find an INS with lying in 1 MK interval. 
Because statistics of the observed sources is poor, we have aver¬ 
aged dP/dT^^ over the age interval from t = 10^ ^ to 10"^ ^ yr (as an 


example). The dashed line is for the HHJ EOS and the mass distri¬ 
bution li (a = 0.1, as in Fig.|^ while the solid line is for the BSk21 
EOS and the mass distribution 2i (a = 0.09, as in Eig.|^. 

Even in the given, rather wide age interval there are only 
six observed sources. They are the neutron star in Pup A (source 
2), PSR J1357-6429 (3), the Vela pulsar (4), PSR B1706-44 (5), 
PSR J0538-f2817 (6), and XMMU J1732-3445 (15). As discussed 
above, sources 2 and 15 can be explained by the cooling theory 
assuming strong proton superfluidity in the neutron star core (e.g., 
[Klochkov et al. 12015} Paper I) which we do not include in our cur¬ 
rent cooling models, for simplicity. Therefore, our present theory is 
relevant only for four sources, 3-6, whose positions are shown in 
Eig.j^by thin vertical dotted lines. With this, extremely poor statis¬ 
tics, we cannot reliably estimate the differential probability distri- 
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log 7;°° 


Figure 8. Theoretical differential probability distributions dPIdT^^ (T^ = 
T^/10^ K) to observe an INS as functions of T^. The probability is av¬ 
eraged over the age interval from t = 10^-^ to 10^-^ yr. The dashed line is 
calculated for the HHJ EOS with the mass distribution li, while the solid 
line is for the BSk21 EOS and the mass distribution 2i. Vertical dotted lines 
show the central surface temperatures for four INSs, sources 3-6. See the 
text for details. 

bution of the observed sources to be directly compared with the¬ 
oretical ones. One can attempt to perform such an analysis using, 
for instance, Bayesian statistics (e.g., |Gelman et al.||2013| l which 
is beyond the scope of this work. Nevertheless, according to Fig. 
I^the obtained theoretical probabilities seem in reasonable accord 
with the data. 

We have plotted a number of other calculated probability dis¬ 
tributions of isolated and accreting neutron stars (for different ages 
t and mass accretion rates (M), for both EOSs, different Md, tr and 
mass distribution functions). We do not present them here. The pos¬ 
sibility of accurate comparison with observations is as questionable 
as in the above example. The statistics of XRTs is better but the val¬ 
ues (M) and extracted from observations are less certain. Our 
theoretical probability distributions seem to agree with the data but 
only qualitatively. Our present qualitative analysis ‘by eye’ seems 
adequate at this stage of the investigation but more advanced anal¬ 
ysis is required in future studies. 

3.4 Less ‘successful’ smaller broadening, 0.05 < cr < 0.08 

Very small broadening of the direct Urea threshold (a smaller than 
0.05) is not acceptable in our model (Section |3.2| ) while the case 

of cr ^ 0.08-0.10 seems most suitable (Section [33] ). Let us 

discuss the intermediate case 0.05 < a < 0.08. In this case, the 
theory can explain the observations but explanations become less 
probable with decreasing a from a ^ 0.08 to cr ^ 0.05 as far as 
mass distributions are concerned as discussed below. On the other 
hand, when a decreases the range of allowable Md widens (Section 

We illustrate these statements in Figs.|^and[^by two exam¬ 
ples for the lowest a ^ 0.05. Fig. [^presents the probability distri¬ 
butions for INSs and XRTs assuming the BSk21 EOS, a = 0.05, 


Md = 1.65 Mg, and the mass functions 3i and 3a (cf. with Figs. 
13 and 1^. The critical mass for opening the direct Urea process is 
almost the same as in our best model for the BSk21 EOS (Fig.[^. 
However, to explain all the data we now have to take narrower mass 
distributions 3i and 3a (Fig.|^ which seem less probable than the 
mass distributions 2i and 2a for the best model. Note also the non¬ 
uniformity of the ‘reference’ curves, especially if compared with 
‘most successful’ BSk21 model (Fig.|^. 

Fig. 1^ shows the probability distributions for INSs and XRTs 
assuming the BSk21 EOS, a = 0.05, and Md = 1.45 Mo. Here, we 
have taken the lowest Md allowed at a = 0.05 (Fig.[^. Again, to 
explain the data we should take rather narrow mass distributions 4i 
and 4a (Fig.|^ shifted to lower masses with respect to the previous 
mass distributions. This model, although formally possible, seems 
rather unlikely, because of three reasons. First, the required mass 
distributions are too narrow; secondly, characteristic neutron star 
masses become uncomfortably low; thirdly, Md is much lower than 
the mass Mdo determined by the BSk21 EOS. Large differences of 
Md from Mdo would be difficult to explain from theoretical point 
of view. 

Therefore, when we decrease a below cr ^ 0.1 but keep Md ~ 
Mdo, we have to narrow neutron star mass distributions (which does 
not seem likely). If we, additionally, take lower Md we increase 
the difference between Md and Mdo and reduce typical masses of 
neutron stars (making our models even less likely). 

Note that the majority of INSs should have masses M < Md. 
On the other hand, some neutron stars in XRTs should be more 
massive than Md (Fig. [^. First of all, this is needed to explain 
the observations of SAX J1804.4-3 658 (source 13;|Campana et ^ 
|2002 1 [Galloway & Cumming|2006[ [Heinke et al.|2009| l which in- 

dicate the operation of direct Urea process. 


3.5 Wider broadening, cr > 0.12 

For completeness of our consideration let us discuss the case of 
relatively large broadening of the direct Urea threshold, a > 0.12. 
The general trend is clear. If we broaden the threshold, the direct 
Urea process operates in wider regions of neutron star cores and 
makes the stars colder. This destroys the agreement of the theory 
with observations. Formally, it is reflected in the fact (Figs. miD 
that at a > 0.12 we obtain MDmin > ^^Dmax- In other words, our 
qualitative analysis of Sectionj^signals that successful solutions of 
the problem disappear. For instance, as seen from figs 15 and 16 of 
Paper I plotted for the HHJ EOS at Md = Mdo and a = 0.2, the 
theory cannot explain the data. 

Finally, we mention that at MDmax <4^Dmin three cases are for¬ 
mally possible, (i) Md < M^max < Momin; (ii) Momax < Mumin < Md 
and (hi) MDmax < Md < MDmin. Case (i) is similar to what is shown 
in Paper I; theoretical domains of ‘warm’ stars (in successful ex¬ 
planations) become cooler but domains of ‘cold’ stars remain un¬ 
changed. Case (ii) is the opposite; ‘warmer’ domains remain un¬ 
changed but ‘colder’ ones become slightly warmer. The last case 
(hi) is a combination of two previous ones; domains of ‘warmer’ 
stars become cooler and ‘colder’ stars warmer. It is unlikely to 
reach satisfactory agreement with the observations of INSs and 
XRTs in all these cases. It would also be difficult to explain rela¬ 
tively large broadening of the direct Urea threshold from theoretical 
point of view. 
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Figure 9. Same as in Figs.|^and[7]but for the BSk21 EOS, a = 0.05, Md = 1.65 Mo, and mass distributions 3i and 3a. 



Figure 10. Same as in Figs.[^[7]and|^but for the BSk21 EOS, a = 0.05, Md = 1.45 Mo, and mass distributions 4i and 4a. 


4 CONCLUSIONS 

We have extended the statistical approach of Paper I to study ther¬ 
mal evolution of middle-aged INSs and old accreting quasistation¬ 
ary neutron stars in XRTs in two respects. 

First, we have considered two EOSs in neutron star cores (the 
HHJ and BSk21 EOSs, Table|^ instead of one (HHJ) EOS in Pa- 
per I. Secondly, we have allowed the theoretical density threshold 
Pd for the direct Urea onset to be not only broadened (as described 
by the parameter a ^ ApD /pdo) but also shifted from its position 
Pdo determined by a given EOS. In this way we have actually ex¬ 
tended a class of studied EOSs because in many cases theoretical 
EOSs give almost the same pressure P{p) and neutron star models 
M{R) but have sufficiently different direct Urea thresholds pd (com¬ 
pare, for instance, the BSkl9, BSk20, and BSk21 EOSs; [Goriely,| 


[Charnel & Pearson|2010[[Pearson et aL|2012[[Potekhin et al.|2013| l 

We have confirmed the principal conclusion of Paper I that 
the theory can explain observations of INSs and XRTs assuming 
the onset of direct Urea process in sufficiently massive stars with 
the broadened density threshold. 

The main conclusions of the present investigation are as fol¬ 
lows. 

(i) At too small (a < 0.05) and too high (a > 0.12) broadening 
a successful explanation of the data is unlikely. 

(ii) At intermediate broadening (0.05 < a < 0.12) one can ex¬ 
plain the data in a certain interval of critical masses Md (MDmin ^ 
Md < MDmax) of stars, where the direct Urea sets in, by tuning the 
mass distributions /i(M) and /a(M) of isolated and accreting neu- 
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tron stars. The results do not depend significantly on EOS (HHJ or 
BSk21). The majority of INSs should have masses M < Mq, but 
some accreting neutron stars should be more massive. 

(iii) The maximum mass Mo^ax ~ 1.7 - 1.8 MG is almost in¬ 
dependent of a, while the minimum mass Momin decreases linearly 
with decreasing a. 

(iv) Most ‘successful’ explanations correspond to cr ^ 0.08 - 
0.10 and Mjy ^ 1.6-1.8 Mo, in which case Momin is close to Momax 
and the mass distributions fi(M) and /a(M) resemble distributions 
1 and 2 in Fig.|^ 

(v) When a decreases from ‘most successful’ values cr ~ 0.1 to 
the lowest values a ^ 0.05, the minimum mass Momin decreases 
and allowable interval of Md becomes wider whereas the mass dis¬ 
tributions fi(M) and /a(M) become uncomfortably narrow (like dis¬ 
tributions 3 and 4 in Fig.|^. If, in addition, Md is taken much below 
^ 1.6 Mo, one gets uncomfortably low neutron star masses. 

(vi) Large shifts of Md from the value Mdo, determined by a 
given EOS, and large broadenings a of the direct Urea threshold, as 
well as too narrow mass distributions /i(M) and /a(M) are unlikely 
from theoretical point of view; such formal solutions have to be 
regarded as less probable. With this in mind, the values a ^ 0.08 - 
0.10 and Md ~ 1.6 —1.8 Mo seem more probable than the others. 

Let us stress that our limitations of Md are not rigorous (rather 
a number of indirect evidences). Our conclusions have to be re¬ 
garded as preliminary and can be refined. First of all, it would be 
instructive to replace our phenomenological model of shifted and 
broadened direct Urea threshold by some physical models taking 
into account the effects of superfluidity and/or magnetic fields as 
well as possible effects of nuclear physics (e.g., |Baiko & Yakovlev| 
|1999| [Yakovlev et al.|[200T] and references therein). In addition, 
one can try different mass distributions of cooling and accreting 
neutron stars (e.g., |Posselt et al.|2008[|Kiziltan et al.|2013| l and for¬ 
malize the comparison of theoretical and observational probability 
distributions of isolated and accreting neutron stars using accurate 
statistical methods (e.g., |Gelman et al.|2013| l instead of comparison 
‘by eye’. It would be good to implement observational selection ef¬ 
fects in this analysis (which seems to be a very complicated prob¬ 
lem). These refinements would hopefully clarify the properties of 
the direct Urea process in neutron stars which are closely related to 
fundamental properties of superdense matter in neutron star cores. 
However, such improvements are beyond the scope of this paper. 
Note that statistical studies of evolution of neutron stars have been 
undertaken previously (e.g., [Popov et al.|2006[|Posselt et al.|2008j ) 
but using quite different approaches. 
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